Modeling catalytic combustion om methane using RMG-Cat¶

Based somewhat on the notebook to accompany the manuscript:
"Automatic generation of microkinetic mechanisms for heterogeneous catalysis" by
C. Franklin Goldsmith, School of Engineering, Brown University, franklin_goldsmith@brown.edu, and
Richard H. West, Department of Chemical Engineering, Northeastern University, r.west@northeastern.edu

As modified to prepare data for the AIChE conference presentation: 304c Incorporation of Linear Scaling Relations into Automatic Mechanism Generation for Heterogeneous Catalysis Richard H. West, Department of Chemical Engineering, Northeastern University, Boston, MA and C. Franklin Goldsmith, Engineering, Brown University, Providence, RI Tuesday, October 31, 2017 https://aiche.confex.com/aiche/2017/meetingapp.cgi/Paper/500172 https://www.slideshare.net/richardhwest/incorporation-of-linear-scaling-relations-into-automatic-mechanism-generation-for-heterogeneous-catalysis

Then further updated to model catalytic combustion of methanol.

First, we print what git commit we were on when we ran this notebook, for both the source code (RMG-Py) and the database.

In [1]:
%%bash
export RMGpy=~/Code/RMG-Py
pwd
git log -n1 --pretty=oneline
cd ../RMG-Py
pwd
git log -n1 --pretty=oneline
cd ../RMG-database
pwd
git log -n1 --pretty=oneline
/Users/emilymazeau/Code/linear-scaling-tests
7eba4f65b6db619f26dec17ee8e6cadef9e24f04 trying again
/Users/emilymazeau/Code/RMG-Py
474969c03bd48274a2af7d1a9b1bbed74ed46e85 Add 'X' to adsorbate names taken from gas phase thermo libraries.
/Users/emilymazeau/Code/RMG-database
48e0a55a7d3f8270c5ffe2dc78f850ed50f7f4f2 Merge remote-tracking branch 'official/master' into cat, again

Model generation¶

We start with a base input file to generate a mechanism for CH3OH plus O2. First we print the input file we'll use to generate the model.

In [2]:
%cat base/input.py
# Data sources
database(
    thermoLibraries=['surfaceThermo', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'],
    reactionLibraries = [('Deutschmann_Ni', False)],
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = 'default',
    kineticsEstimator = 'rate rules',
    bindingEnergies = { # default values for Ni(111)
                       'C':(-5.997, 'eV/molecule'),
                       'H':(-2.778, 'eV/molecule'),
                       'O':(-4.485, 'eV/molecule'),
                       }
)

# List of species
species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)
species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='N2',
    reactive=False,
    structure=SMILES("N#N"),
)

species(
    label='CO2',
    reactive=True,
    structure=SMILES("O=C=O"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

species(
    label='H2',
    reactive=True,
    structure=SMILES("[H][H]"),
)

species(
    label='CO',
    reactive=True,
    structure=SMILES("[C-]#[O+]"),
)

species(
    label='C2H6',
    reactive=True,
    structure=SMILES("CC"),
)

species(
    label='CH2O',
    reactive=True,
    structure=SMILES("C=O"),
)

species(
    label='CH3',
    reactive=True,
    structure=SMILES("[CH3]"),
)

species(
    label='C3H8',
    reactive=True,
    structure=SMILES("CCC"),
)

species(
    label='H',
    reactive=True,
    structure=SMILES("[H]"),
)

species(
    label='C2H5',
    reactive=True,
    structure=SMILES("C[CH2]"),
)

species(
    label='CH3OH',
    reactive=True,
    structure=SMILES("CO"),
)

species(
    label='HCO',
    reactive=True,
    structure=SMILES("[CH]=O"),
)

species(
    label='CH3CHO',
    reactive=True,
    structure=SMILES("CC=O"),
)

species(
    label='OH',
    reactive=True,
    structure=SMILES("[OH]"),
)

species(
    label='C2H4',
    reactive=True,
    structure=SMILES("C=C"),
)

#----------
# Reaction systems
surfaceReactor(
    temperature=(1200,'K'),
    initialPressure=(1.01325, 'bar'),
    initialGasMoleFractions={
        "CH3OH": 0.11764706,
        "O2": 0.17647059,
        "N2": 0.70588235,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    surfaceSiteDensity=(2.9e-9, 'mol/cm^2'),
#    terminationConversion = { "CH4":0.9,},
    terminationTime=(1.0, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=1e-2,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000
)

options(
    units='si',
    saveRestartPeriod=None,
    generateOutputHTML=True,
    generatePlots=False,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

Then we try running it. This will take Richard's computer about four minutes and Emily's computer 25 minutes.

In [3]:
%%bash
python ~/Code/RMG-Py/rmg.py base/input.py > /dev/null
tail -n12 base/RMG.log
    Execution time (DD:HH:MM:SS): 00:00:23:43
    Memory used: 1050.95 MB

Making seed mechanism...
Performing final model checks...

MODEL GENERATION COMPLETED

The final model core has 101 species and 324 reactions
The final model edge has 1088 species and 2975 reactions

RMG execution terminated at Mon Jul  9 15:32:04 2018
:root:Removing old /Users/emilymazeau/Code/linear-scaling-tests/base/RMG_backup.log
:root:Moving /Users/emilymazeau/Code/linear-scaling-tests/base/RMG.log to /Users/emilymazeau/Code/linear-scaling-tests/base/RMG_backup.log

/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py:1135: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  data = item.generateReverseRateCoefficient()
/Users/emilymazeau/Code/RMG-Py/rmgpy/thermo/thermoengine.py:58: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  wilhoit = thermo0.toWilhoit(B=1000.)
/Users/emilymazeau/.local/lib/python2.7/site-packages/scipy/optimize/optimize.py:1742: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  fx = func(x, *args)
/Users/emilymazeau/.local/lib/python2.7/site-packages/scipy/optimize/optimize.py:1794: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  fu = func(x, *args)
/Users/emilymazeau/.local/lib/python2.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

There are 70 species and 218 reactions (?)

Data processing¶

Next we will import some libraries and set things up to start importing and analyzing the simulation results.

In [4]:
%config InlineBackend.figure_format = 'retina'
In [5]:
%matplotlib inline
from matplotlib import pyplot as plt
import matplotlib

# The default output Type 3 (Type3) fonts can't be edited in Adobe Illustrator
# but Type 42 (TrueType) fonts can be, making it easier to move labels around
# slightly to improve layout before publication.
matplotlib.rcParams['pdf.fonttype'] = 42 

# Seaborn helps make matplotlib graphics nicer
import seaborn as sns
sns.set_style("ticks")
sns.set_context("paper", font_scale=1.5, rc={"lines.linewidth": 2.0})

import os
import re
import pandas as pd
import numpy as np
import shutil
import subprocess
import multiprocessing
In [6]:
def get_last_csv_file(job_directory):
    """
    Find the CSV file from the largest model in the provided job directory.
    
    For CSV files named `simulation_1_13.csv` you want 13 to be the highest number.
    """
    solver_directory = os.path.join(job_directory,'solver')
    csv_files = sorted([f for f in os.listdir(solver_directory) if f.endswith('.csv') ],
                       key=lambda f: int(f[:-4].split('_')[2]))
    return os.path.join(solver_directory, csv_files[-1])
    
job_directory = 'base'
get_last_csv_file(job_directory)
Out[6]:
'base/solver/simulation_1_101.csv'

We will use Pandas to import the csv file

In [7]:
last_csv_file = get_last_csv_file(job_directory)
data = pd.read_csv(last_csv_file)
data
Out[7]:
Time (s) Volume (m^3) N2 Ar He Ne X(1) CH4(2) O2(3) CO2(4) ... SX(943) SX(974) C2H3X(50) SX(617) C3H8X(43) C3H7X(41) SX(254) C3H7X(42) SX(169) C3H7(62)
0 0.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.285560 0.000000e+00 1.764706e-01 0.000000e+00 ... 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1 1.856202e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.717598e-50 1.764706e-01 1.680474e-101 ... 1.815842e-157 -9.103578e-194 2.978876e-81 5.589589e-89 -8.244069e-106 -7.465910e-108 -1.239111e-158 -2.488637e-108 -1.757593e-143 -1.956214e-116
2 5.197367e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.129194e-49 1.764706e-01 1.856996e-99 ... 5.575303e-155 1.005338e-190 6.663628e-80 2.893332e-87 1.383778e-104 1.253163e-106 3.173318e-156 4.177209e-107 3.571228e-141 3.283537e-115
3 1.187970e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.669248e-48 1.764706e-01 1.680071e-97 ... 5.325582e-152 4.394961e-187 1.316621e-78 1.263342e-85 1.631208e-102 1.477238e-104 3.592568e-153 4.924127e-105 1.980076e-138 3.870682e-113
4 2.524435e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.049355e-47 1.764706e-01 1.271504e-95 ... 3.638884e-149 1.276425e-183 2.342857e-77 4.703106e-84 1.289394e-100 1.167688e-102 2.704337e-150 3.892294e-103 7.259774e-136 3.059627e-111
5 5.197367e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.030870e-46 1.764706e-01 8.833452e-94 ... 2.142735e-146 3.095097e-180 3.951573e-76 1.620291e-82 9.125113e-99 8.263791e-101 1.667862e-147 2.754597e-101 2.208166e-133 2.165361e-109
6 1.054323e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.986439e-45 1.764706e-01 5.887460e-92 ... 1.175272e-143 6.887802e-177 6.490479e-75 5.377156e-81 6.135928e-97 5.556756e-99 9.355449e-145 1.852252e-99 6.150769e-131 1.456100e-107
7 2.123496e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.054889e-44 1.764706e-01 3.844737e-90 ... 6.226505e-141 1.469950e-173 1.052137e-73 1.752068e-79 4.024417e-95 3.644551e-97 5.011346e-142 1.214850e-97 1.641749e-128 9.551058e-106
8 4.261841e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.270415e-43 1.764706e-01 2.485498e-88 ... 3.242679e-138 3.072856e-170 1.694442e-72 5.657321e-78 2.607249e-93 2.361150e-95 2.624126e-139 7.870500e-96 4.291093e-126 6.188796e-104
9 8.538531e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.626975e-42 1.764706e-01 1.598725e-86 ... 1.674409e-135 6.357963e-167 2.719981e-71 1.818495e-76 1.678830e-91 1.520365e-93 1.358697e-136 5.067884e-94 1.109958e-123 3.986399e-102
10 1.709191e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.105846e-41 1.764706e-01 1.025754e-84 ... 8.609417e-133 1.308786e-163 4.359146e-70 5.832261e-75 1.077724e-89 9.759975e-92 6.995577e-134 3.253325e-92 2.856230e-121 2.560837e-100
11 3.419867e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.686384e-40 1.764706e-01 6.573068e-83 ... 4.417375e-130 2.687252e-160 6.980580e-69 1.868418e-73 6.907926e-88 6.255885e-90 3.591769e-131 2.085295e-90 7.330871e-119 1.643698e-98
12 6.841220e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.349791e-39 1.764706e-01 4.209408e-81 ... 2.264094e-127 5.510526e-157 1.117433e-67 5.982292e-72 4.424433e-86 4.006809e-88 1.841561e-128 1.335603e-88 1.879129e-116 1.055674e-96
13 1.368392e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.080106e-38 1.764706e-01 2.694878e-79 ... 1.159831e-124 1.129277e-153 1.788530e-66 1.914870e-70 2.832713e-84 2.565332e-86 9.435417e-126 8.551107e-87 4.813678e-114 6.796104e-95
14 2.736933e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.641944e-38 1.764706e-01 1.725013e-77 ... 5.939906e-122 2.313497e-150 2.862815e-65 6.128444e-69 1.813280e-82 1.642124e-84 4.832665e-123 5.473748e-85 1.232699e-111 4.397993e-93
15 5.474015e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.913992e-37 1.764706e-01 1.104137e-75 ... 3.041635e-119 4.738796e-147 4.583539e-64 1.961239e-67 1.160610e-80 1.051059e-82 2.474824e-120 3.503531e-83 3.156219e-109 2.875993e-91
16 1.094818e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.531368e-36 1.764706e-01 7.067676e-74 ... 1.557420e-116 9.705825e-144 7.342818e-63 6.276186e-66 7.428253e-79 6.727098e-81 1.267334e-117 2.242366e-81 8.080571e-107 1.918820e-89
17 2.189651e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.425164e-35 1.764706e-01 4.525567e-72 ... 7.974261e-114 1.987830e-140 1.177736e-61 2.008415e-64 4.754193e-77 4.305444e-79 6.491131e-115 1.435148e-79 2.068708e-104 1.328041e-87
18 4.379316e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.540157e-34 1.764706e-01 2.901684e-70 ... 4.082893e-111 4.071150e-137 1.893572e-60 6.426984e-63 3.042718e-75 2.755515e-77 3.327583e-112 9.185051e-78 5.295991e-102 9.779147e-86
19 8.758647e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.832135e-33 1.764706e-01 1.870304e-68 ... 2.090462e-108 8.337765e-134 3.059113e-59 2.056644e-61 1.947349e-73 1.763539e-75 1.711788e-109 5.878462e-76 1.355783e-99 7.896569e-84
20 1.751731e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.265708e-32 1.764706e-01 1.230383e-66 ... 1.070325e-105 1.707569e-130 4.988628e-58 6.581274e-60 1.246305e-71 1.128666e-73 2.758656e-106 3.762219e-74 7.502154e-97 7.150303e-82
21 3.503463e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.812562e-31 1.764706e-01 8.721977e-65 ... 5.480112e-103 3.497053e-127 8.282737e-57 2.106010e-58 7.976329e-70 7.223442e-72 1.565004e-103 2.407814e-72 1.924569e-94 7.259688e-80
22 7.006928e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.450040e-30 1.764706e-01 7.735366e-63 ... 2.805860e-100 7.161743e-124 1.421534e-55 6.739236e-57 5.104812e-68 4.622969e-70 1.107687e-100 1.540989e-70 4.926853e-92 8.081027e-78
23 1.051039e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.304789e-30 1.764706e-01 4.904492e-62 ... 2.412932e-99 7.587052e-123 5.421762e-55 3.143839e-56 2.914548e-67 2.639443e-69 6.051405e-100 8.798144e-70 2.678597e-91 5.270121e-77
24 1.751732e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.016426e-29 1.764706e-01 2.662649e-60 ... 4.863211e-97 6.210470e-120 4.684244e-54 4.212949e-55 7.541711e-66 6.829850e-68 4.802546e-97 2.276616e-68 4.070377e-89 2.238711e-75
25 2.382356e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.653069e-29 1.764706e-01 1.413330e-59 ... 3.609156e-96 6.148633e-119 1.358069e-53 1.398248e-54 3.129736e-65 2.834321e-67 5.560508e-96 9.447736e-68 2.904789e-88 1.025904e-74
26 3.012980e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.959133e-29 1.764706e-01 5.394133e-59 ... 1.516579e-95 3.060522e-118 3.232480e-53 3.757111e-54 9.856511e-65 8.926157e-67 2.028853e-95 2.975385e-67 9.721032e-88 3.701024e-74
27 3.643603e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.523261e-28 1.764706e-01 1.655002e-58 ... 5.483267e-95 1.323297e-117 6.524392e-53 8.311626e-54 2.518899e-64 2.281141e-66 1.072338e-94 7.603801e-67 3.503235e-87 1.042595e-73
28 4.904851e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.367896e-28 1.764706e-01 9.821801e-58 ... 5.328879e-94 2.029880e-116 1.820282e-52 2.615099e-53 1.026388e-63 9.295076e-66 2.191740e-93 3.098358e-66 3.159397e-86 4.781406e-73
29 6.166098e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.309620e-28 1.764706e-01 3.769097e-57 ... 2.428248e-93 1.166697e-115 4.120776e-52 6.468490e-53 3.077094e-63 2.786648e-65 1.239401e-92 9.288826e-66 1.327957e-85 1.575102e-72
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
18095 9.448643e-01 0.098469 0.705882 0.0 0.0 0.0 0.202383 2.648879e-10 2.247721e-12 8.470734e-02 ... 8.884909e-08 3.756688e-17 3.320348e-06 5.720220e-04 2.445150e-11 6.351454e-14 3.104531e-11 1.318308e-13 8.950091e-17 1.223302e-10
18096 9.468529e-01 0.098469 0.705882 0.0 0.0 0.0 0.202405 2.632985e-10 2.245945e-12 8.471634e-02 ... 8.881663e-08 3.753832e-17 3.321653e-06 5.720207e-04 2.445413e-11 6.350703e-14 3.100971e-11 1.318152e-13 8.934495e-17 1.223025e-10
18097 9.488415e-01 0.098469 0.705882 0.0 0.0 0.0 0.202427 2.617215e-10 2.244176e-12 8.472532e-02 ... 8.878429e-08 3.750987e-17 3.322954e-06 5.720194e-04 2.445675e-11 6.349955e-14 3.097414e-11 1.317997e-13 8.918996e-17 1.222750e-10
18098 9.508301e-01 0.098469 0.705882 0.0 0.0 0.0 0.202448 2.601567e-10 2.242413e-12 8.473427e-02 ... 8.875206e-08 3.748153e-17 3.324252e-06 5.720181e-04 2.445937e-11 6.349209e-14 3.093862e-11 1.317842e-13 8.903595e-17 1.222475e-10
18099 9.528186e-01 0.098469 0.705882 0.0 0.0 0.0 0.202470 2.586040e-10 2.240658e-12 8.474320e-02 ... 8.871995e-08 3.745330e-17 3.325547e-06 5.720168e-04 2.446197e-11 6.348465e-14 3.090313e-11 1.317687e-13 8.888289e-17 1.222202e-10
18100 9.548072e-01 0.098469 0.705882 0.0 0.0 0.0 0.202492 2.570635e-10 2.238909e-12 8.475209e-02 ... 8.868796e-08 3.742518e-17 3.326838e-06 5.720155e-04 2.446457e-11 6.347723e-14 3.086769e-11 1.317533e-13 8.873079e-17 1.221929e-10
18101 9.567958e-01 0.098469 0.705882 0.0 0.0 0.0 0.202513 2.555349e-10 2.237167e-12 8.476096e-02 ... 8.865608e-08 3.739717e-17 3.328125e-06 5.720142e-04 2.446716e-11 6.346984e-14 3.083229e-11 1.317380e-13 8.857964e-17 1.221657e-10
18102 9.587844e-01 0.098469 0.705882 0.0 0.0 0.0 0.202535 2.540181e-10 2.235431e-12 8.476980e-02 ... 8.862432e-08 3.736927e-17 3.329410e-06 5.720129e-04 2.446975e-11 6.346246e-14 3.079692e-11 1.317227e-13 8.842943e-17 1.221386e-10
18103 9.607729e-01 0.098469 0.705882 0.0 0.0 0.0 0.202556 2.525132e-10 2.233702e-12 8.477861e-02 ... 8.859266e-08 3.734147e-17 3.330691e-06 5.720117e-04 2.447232e-11 6.345511e-14 3.076160e-11 1.317074e-13 8.828015e-17 1.221115e-10
18104 9.627615e-01 0.098469 0.705882 0.0 0.0 0.0 0.202577 2.510199e-10 2.231980e-12 8.478739e-02 ... 8.856112e-08 3.731378e-17 3.331968e-06 5.720104e-04 2.447489e-11 6.344778e-14 3.072631e-11 1.316922e-13 8.813180e-17 1.220846e-10
18105 9.647501e-01 0.098469 0.705882 0.0 0.0 0.0 0.202598 2.495382e-10 2.230264e-12 8.479615e-02 ... 8.852969e-08 3.728620e-17 3.333242e-06 5.720091e-04 2.447745e-11 6.344047e-14 3.069107e-11 1.316770e-13 8.798436e-17 1.220577e-10
18106 9.667387e-01 0.098469 0.705882 0.0 0.0 0.0 0.202620 2.480679e-10 2.228555e-12 8.480488e-02 ... 8.849837e-08 3.725872e-17 3.334513e-06 5.720078e-04 2.448001e-11 6.343318e-14 3.065587e-11 1.316619e-13 8.783784e-17 1.220309e-10
18107 9.687273e-01 0.098469 0.705882 0.0 0.0 0.0 0.202641 2.466091e-10 2.226852e-12 8.481358e-02 ... 8.846715e-08 3.723134e-17 3.335781e-06 5.720066e-04 2.448255e-11 6.342591e-14 3.062070e-11 1.316468e-13 8.769221e-17 1.220042e-10
18108 9.707158e-01 0.098469 0.705882 0.0 0.0 0.0 0.202662 2.451615e-10 2.225155e-12 8.482226e-02 ... 8.843605e-08 3.720407e-17 3.337045e-06 5.720053e-04 2.448509e-11 6.341866e-14 3.058558e-11 1.316318e-13 8.754749e-17 1.219776e-10
18109 9.727044e-01 0.098469 0.705882 0.0 0.0 0.0 0.202683 2.437251e-10 2.223464e-12 8.483091e-02 ... 8.840505e-08 3.717690e-17 3.338306e-06 5.720040e-04 2.448762e-11 6.341143e-14 3.055049e-11 1.316168e-13 8.740365e-17 1.219511e-10
18110 9.746930e-01 0.098469 0.705882 0.0 0.0 0.0 0.202704 2.422999e-10 2.221780e-12 8.483953e-02 ... 8.837416e-08 3.714983e-17 3.339564e-06 5.720028e-04 2.449015e-11 6.340422e-14 3.051545e-11 1.316018e-13 8.726069e-17 1.219246e-10
18111 9.766816e-01 0.098469 0.705882 0.0 0.0 0.0 0.202725 2.408857e-10 2.220102e-12 8.484813e-02 ... 8.834338e-08 3.712286e-17 3.340818e-06 5.720015e-04 2.449266e-11 6.339704e-14 3.048044e-11 1.315869e-13 8.711861e-17 1.218982e-10
18112 9.786701e-01 0.098469 0.705882 0.0 0.0 0.0 0.202745 2.394824e-10 2.218430e-12 8.485670e-02 ... 8.831270e-08 3.709600e-17 3.342070e-06 5.720003e-04 2.449517e-11 6.338987e-14 3.044548e-11 1.315720e-13 8.697740e-17 1.218719e-10
18113 9.806587e-01 0.098469 0.705882 0.0 0.0 0.0 0.202766 2.380899e-10 2.216764e-12 8.486525e-02 ... 8.828212e-08 3.706923e-17 3.343318e-06 5.719990e-04 2.449768e-11 6.338272e-14 3.041055e-11 1.315572e-13 8.683705e-17 1.218457e-10
18114 9.826473e-01 0.098469 0.705882 0.0 0.0 0.0 0.202787 2.367082e-10 2.215104e-12 8.487377e-02 ... 8.825165e-08 3.704256e-17 3.344563e-06 5.719978e-04 2.450017e-11 6.337559e-14 3.037566e-11 1.315424e-13 8.669755e-17 1.218196e-10
18115 9.846359e-01 0.098469 0.705882 0.0 0.0 0.0 0.202808 2.353371e-10 2.213450e-12 8.488227e-02 ... 8.822128e-08 3.701598e-17 3.345805e-06 5.719966e-04 2.450266e-11 6.336849e-14 3.034082e-11 1.315276e-13 8.655890e-17 1.217935e-10
18116 9.866244e-01 0.098469 0.705882 0.0 0.0 0.0 0.202828 2.339767e-10 2.211803e-12 8.489074e-02 ... 8.819101e-08 3.698951e-17 3.347044e-06 5.719953e-04 2.450515e-11 6.336140e-14 3.030601e-11 1.315129e-13 8.642110e-17 1.217675e-10
18117 9.886130e-01 0.098469 0.705882 0.0 0.0 0.0 0.202849 2.326267e-10 2.210161e-12 8.489918e-02 ... 8.816084e-08 3.696313e-17 3.348280e-06 5.719941e-04 2.450762e-11 6.335433e-14 3.027124e-11 1.314982e-13 8.628412e-17 1.217416e-10
18118 9.906016e-01 0.098469 0.705882 0.0 0.0 0.0 0.202869 2.312871e-10 2.208524e-12 8.490761e-02 ... 8.813078e-08 3.693684e-17 3.349513e-06 5.719928e-04 2.451009e-11 6.334728e-14 3.023651e-11 1.314836e-13 8.614798e-17 1.217158e-10
18119 9.925902e-01 0.098469 0.705882 0.0 0.0 0.0 0.202890 2.299578e-10 2.206894e-12 8.491600e-02 ... 8.810081e-08 3.691065e-17 3.350743e-06 5.719916e-04 2.451256e-11 6.334025e-14 3.020182e-11 1.314690e-13 8.601265e-17 1.216900e-10
18120 9.945788e-01 0.098469 0.705882 0.0 0.0 0.0 0.202910 2.286388e-10 2.205270e-12 8.492438e-02 ... 8.807094e-08 3.688455e-17 3.351969e-06 5.719904e-04 2.451501e-11 6.333324e-14 3.016717e-11 1.314545e-13 8.587814e-17 1.216643e-10
18121 9.965673e-01 0.098469 0.705882 0.0 0.0 0.0 0.202930 2.273299e-10 2.203651e-12 8.493272e-02 ... 8.804117e-08 3.685855e-17 3.353193e-06 5.719892e-04 2.451746e-11 6.332625e-14 3.013256e-11 1.314400e-13 8.574444e-17 1.216387e-10
18122 9.985559e-01 0.098469 0.705882 0.0 0.0 0.0 0.202951 2.260311e-10 2.202038e-12 8.494105e-02 ... 8.801150e-08 3.683264e-17 3.354414e-06 5.719879e-04 2.451990e-11 6.331927e-14 3.009798e-11 1.314255e-13 8.561154e-17 1.216131e-10
18123 1.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.202965 2.250942e-10 2.200870e-12 8.494708e-02 ... 8.799001e-08 3.681388e-17 3.355299e-06 5.719871e-04 2.452167e-11 6.331422e-14 3.007290e-11 1.314150e-13 8.551553e-17 1.215946e-10
18124 1.000544e+00 0.098469 0.705882 0.0 0.0 0.0 0.202971 2.247423e-10 2.200430e-12 8.494935e-02 ... 8.798192e-08 3.680682e-17 3.355632e-06 5.719867e-04 2.452234e-11 6.331232e-14 3.006345e-11 1.314110e-13 8.547943e-17 1.215877e-10

18125 rows × 103 columns

In [8]:
def get_pandas_data(job_directory):
    """
    Get the last CSV file from the provided job directory,
    import it into a Pandas data frame, and tidy it up a bit.
    """
    last_csv_file = get_last_csv_file(job_directory)
    data = pd.read_csv(last_csv_file)
    
    # Make the Time into an index and remove it as a column
    data.index = data['Time (s)']
    del data['Time (s)']
    # Remove inerts that RMG added automatically but we're not using
    for i in 'Ar He Ne'.split():
        del data[i]
    # Remove the Volume column
    del data['Volume (m^3)']
    # Set any amounts that are slightly negative (within the ODE solver's ATOL tolerance), equal to zero
    # to allow 'area' plots without warnings or errors.
    # Anything more negative than -1e-16 probably indicates a bug in RMG and should not be hidden here ;-)
    data[(data<0) & (data>-1e-16)] = 0
    return data
In [9]:
def rename_columns(data):
    """
    Removes the number (##) from the end of the column names, in place,
    unless doing so would make the names collide.
    Also renames a few species so the plot labels match the names in the manuscript.
    """
    import re
    old = data.columns
    new = [re.sub('\(\d+\)$','',n) for n in old]
    # don't translate names that would no longer be unique
    mapping = {k:v for k,v in zip(old,new) if new.count(v)==1}
    data.rename(columns=mapping, inplace=True)
    
    # Now change a few species that are named differently in the manuscript
    # than in the thermodynamics database used to build the model,
    # so that the plot labels match the manuscript.
    mapping = {'COX':'COvdwX', 'OCX': 'CO=X', 'C2H3X':'CH3CX', 'C2H3OX':'CH3CXO'}
    data.rename(columns=mapping, inplace=True)
In [10]:
data1 = get_pandas_data('base')
rename_columns(data1)
data1.columns
Out[10]:
Index([u'N2', u'X', u'CH4', u'O2', u'CO2', u'H2O', u'H2', u'CO', u'C2H6',
       u'CH2O', u'CH3', u'C3H8', u'H', u'C2H5', u'CH3OH', u'HCO', u'CH3CHO',
       u'OH', u'C2H4(18)', u'CH3OHX', u'OX', u'HX', u'CH3OX(45)', u'OCCO',
       u'HOX', u'CH2OX(135)', u'CHOX', u'H2X', u'H2OX', u'HO2X', u'SX(156)',
       u'SX(155)', u'CX', u'CO=X', u'COvdwX', u'CH3OX(44)', u'SX(179)',
       u'CH3X', u'CH2X', u'CHX', u'CXHO', u'SX(165)', u'CH4X', u'CO2X',
       u'SX(185)', u'SX(183)', u'SX(246)', u'SX(225)', u'SX(143)', u'CH2OH',
       u'HO2', u'CH2OX(40)', u'HOCXO', u'SX(265)', u'CH3O', u'SX(191)',
       u'SX(333)', u'SX(356)', u'SX(178)', u'HOOH', u'HOOHX', u'CH2(S)',
       u'C2H6O', u'SX(422)', u'C2H4X', u'C2H4(32)', u'CHO4', u'O(T)', u'CH3O3',
       u'DME', u'DMEX', u'C3H6O', u'S(439)', u'C2H4O(412)', u'C2H4O(30)',
       u'C2H4OX', u'C2H4O(139)', u'SX(620)', u'S(411)', u'S(692)', u'S(693)',
       u'C2H6X', u'CH4O2', u'SX(878)', u'SX(171)', u'SX(885)', u'SX(914)',
       u'SX(161)', u'SX(943)', u'SX(974)', u'CH3CX', u'SX(617)', u'C3H8X',
       u'C3H7X(41)', u'SX(254)', u'C3H7X(42)', u'SX(169)', u'C3H7'],
      dtype='object')
In [11]:
# Test it with some plots
data1[['CH3OH', 'O2']].plot.line()
data1[['CO', 'H2']].plot.line()
data1[['H2O']].plot.line()
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x104a0ea90>
In [12]:
species_names = data1.columns
# just the gas phase species that aren't always zero:
gas_phase = [n for n in species_names if 'X' not in n and (data1[n]>0).any()]
# all the surface species
surface_phase = [n for n in species_names if 'X' in n]
surface_phase.remove('X')
data1[gas_phase].plot.line()
data1[surface_phase].plot.line()
Out[12]:
<matplotlib.axes._subplots.AxesSubplot at 0x10a786f50>
In [13]:
print "Significant species (those that exceed 0.001 mol at some point)"
significant = [n for n in data1.columns if(data1[n]>0.001).any()]
with sns.color_palette("hls", len(significant)):
    data1[significant].plot.area(legend='reverse')
Significant species (those that exceed 0.001 mol at some point)
In [14]:
lim = 1e-4
surface = [n for n in data1.columns if 'X' in n and n!='X' and (data1[n]>lim).any() ]
print "The {} surface species that exceed {:.1e} mol at some point".format(len(surface), lim)
total_sites = max(data1['X'])
with sns.color_palette('Set3',len(surface)):
    (data1[surface]/total_sites).plot.area(legend='reverse')
    plt.ylabel('site fraction')
    plt.xlim(0,1e-2) ## notice this is much less than 1 seconde
The 14 surface species that exceed 1.0e-04 mol at some point

Effect of binding energies¶

Now we will use that template 'base' input file to create a ton of other input files with different binding energies, then run them all in RMG-Cat, then process the results.

In [15]:
# First, make a series of input files in separate directories

with open(os.path.join('base', 'input.py')) as infile:
    input_file = infile.read()
    
base_directory = 'binding_energies'
def directory(carbon,oxygen):
    return os.path.join(base_directory, "c{:.3f}o{:.3f}".format(carbon,oxygen))

def make_input(binding_energies):
    """
    Make an input file for the given (carbon,oxygen) tuple (or iterable) of binding energies
    and return the name of the directory in which it is saved.
    """
    carbon, oxygen = binding_energies
    output = input_file
    out_dir = directory(carbon, oxygen)
    carbon_string = "'C':({:f}, 'eV/molecule')".format(carbon)
    output = re.sub("'C':\(.*?, 'eV/molecule'\)", carbon_string, output)
    oxygen_string = "'O':({:f}, 'eV/molecule')".format(oxygen)
    output = re.sub("'O':\(.*?, 'eV/molecule'\)", oxygen_string, output)
    os.path.exists(out_dir) or os.makedirs(out_dir)
    out_file = os.path.join(out_dir, 'input.py')
    with open(out_file,'w') as outfile:
        outfile.write(output)
    shutil.copy(os.path.join('base','run.sh'), out_dir)
    return out_dir

    
print make_input((-8,-3.5))
    
binding_energies/c-8.000o-3.500
In [18]:
def run_simulation(binding_energies):
    """
    Assuming a job file already exists, run it.  This one is local.
    Takes a tuple of binding energies, (carbon, oxygen)
    """
    carbon, oxygen = binding_energies
    job_directory = directory(carbon, oxygen)
    print job_directory
    assert os.path.exists(job_directory)
    return subprocess.check_call('./run.sh', cwd=job_directory)
# a test experiment = (-8,-3.5) make_input(experiment) run_simulation(experiment)
In [19]:
# Revised range
carbon_range = (-7.5, -2)
oxygen_range = (-6.5, -1.5)
grid_size = 9
mesh  = np.mgrid[carbon_range[0]:carbon_range[1]:grid_size*1j, 
                 oxygen_range[0]:oxygen_range[1]:grid_size*1j]

with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
plt.show()
    
mesh
Out[19]:
array([[[-7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   ,
         -7.5   , -7.5   ],
        [-6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125,
         -6.8125, -6.8125],
        [-6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 ,
         -6.125 , -6.125 ],
        [-5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375,
         -5.4375, -5.4375],
        [-4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  ,
         -4.75  , -4.75  ],
        [-4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625,
         -4.0625, -4.0625],
        [-3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 ,
         -3.375 , -3.375 ],
        [-2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875,
         -2.6875, -2.6875],
        [-2.    , -2.    , -2.    , -2.    , -2.    , -2.    , -2.    ,
         -2.    , -2.    ]],

       [[-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ]]])
In [20]:
experiments = mesh.reshape((2,-1)).T
experiments
Out[20]:
array([[-7.5   , -6.5   ],
       [-7.5   , -5.875 ],
       [-7.5   , -5.25  ],
       [-7.5   , -4.625 ],
       [-7.5   , -4.    ],
       [-7.5   , -3.375 ],
       [-7.5   , -2.75  ],
       [-7.5   , -2.125 ],
       [-7.5   , -1.5   ],
       [-6.8125, -6.5   ],
       [-6.8125, -5.875 ],
       [-6.8125, -5.25  ],
       [-6.8125, -4.625 ],
       [-6.8125, -4.    ],
       [-6.8125, -3.375 ],
       [-6.8125, -2.75  ],
       [-6.8125, -2.125 ],
       [-6.8125, -1.5   ],
       [-6.125 , -6.5   ],
       [-6.125 , -5.875 ],
       [-6.125 , -5.25  ],
       [-6.125 , -4.625 ],
       [-6.125 , -4.    ],
       [-6.125 , -3.375 ],
       [-6.125 , -2.75  ],
       [-6.125 , -2.125 ],
       [-6.125 , -1.5   ],
       [-5.4375, -6.5   ],
       [-5.4375, -5.875 ],
       [-5.4375, -5.25  ],
       [-5.4375, -4.625 ],
       [-5.4375, -4.    ],
       [-5.4375, -3.375 ],
       [-5.4375, -2.75  ],
       [-5.4375, -2.125 ],
       [-5.4375, -1.5   ],
       [-4.75  , -6.5   ],
       [-4.75  , -5.875 ],
       [-4.75  , -5.25  ],
       [-4.75  , -4.625 ],
       [-4.75  , -4.    ],
       [-4.75  , -3.375 ],
       [-4.75  , -2.75  ],
       [-4.75  , -2.125 ],
       [-4.75  , -1.5   ],
       [-4.0625, -6.5   ],
       [-4.0625, -5.875 ],
       [-4.0625, -5.25  ],
       [-4.0625, -4.625 ],
       [-4.0625, -4.    ],
       [-4.0625, -3.375 ],
       [-4.0625, -2.75  ],
       [-4.0625, -2.125 ],
       [-4.0625, -1.5   ],
       [-3.375 , -6.5   ],
       [-3.375 , -5.875 ],
       [-3.375 , -5.25  ],
       [-3.375 , -4.625 ],
       [-3.375 , -4.    ],
       [-3.375 , -3.375 ],
       [-3.375 , -2.75  ],
       [-3.375 , -2.125 ],
       [-3.375 , -1.5   ],
       [-2.6875, -6.5   ],
       [-2.6875, -5.875 ],
       [-2.6875, -5.25  ],
       [-2.6875, -4.625 ],
       [-2.6875, -4.    ],
       [-2.6875, -3.375 ],
       [-2.6875, -2.75  ],
       [-2.6875, -2.125 ],
       [-2.6875, -1.5   ],
       [-2.    , -6.5   ],
       [-2.    , -5.875 ],
       [-2.    , -5.25  ],
       [-2.    , -4.625 ],
       [-2.    , -4.    ],
       [-2.    , -3.375 ],
       [-2.    , -2.75  ],
       [-2.    , -2.125 ],
       [-2.    , -1.5   ]])
In [21]:
with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
    plt.plot(*experiments.T, marker='o', linestyle='none')
In [22]:
map(make_input, experiments)
Out[22]:
['binding_energies/c-7.500o-6.500',
 'binding_energies/c-7.500o-5.875',
 'binding_energies/c-7.500o-5.250',
 'binding_energies/c-7.500o-4.625',
 'binding_energies/c-7.500o-4.000',
 'binding_energies/c-7.500o-3.375',
 'binding_energies/c-7.500o-2.750',
 'binding_energies/c-7.500o-2.125',
 'binding_energies/c-7.500o-1.500',
 'binding_energies/c-6.812o-6.500',
 'binding_energies/c-6.812o-5.875',
 'binding_energies/c-6.812o-5.250',
 'binding_energies/c-6.812o-4.625',
 'binding_energies/c-6.812o-4.000',
 'binding_energies/c-6.812o-3.375',
 'binding_energies/c-6.812o-2.750',
 'binding_energies/c-6.812o-2.125',
 'binding_energies/c-6.812o-1.500',
 'binding_energies/c-6.125o-6.500',
 'binding_energies/c-6.125o-5.875',
 'binding_energies/c-6.125o-5.250',
 'binding_energies/c-6.125o-4.625',
 'binding_energies/c-6.125o-4.000',
 'binding_energies/c-6.125o-3.375',
 'binding_energies/c-6.125o-2.750',
 'binding_energies/c-6.125o-2.125',
 'binding_energies/c-6.125o-1.500',
 'binding_energies/c-5.438o-6.500',
 'binding_energies/c-5.438o-5.875',
 'binding_energies/c-5.438o-5.250',
 'binding_energies/c-5.438o-4.625',
 'binding_energies/c-5.438o-4.000',
 'binding_energies/c-5.438o-3.375',
 'binding_energies/c-5.438o-2.750',
 'binding_energies/c-5.438o-2.125',
 'binding_energies/c-5.438o-1.500',
 'binding_energies/c-4.750o-6.500',
 'binding_energies/c-4.750o-5.875',
 'binding_energies/c-4.750o-5.250',
 'binding_energies/c-4.750o-4.625',
 'binding_energies/c-4.750o-4.000',
 'binding_energies/c-4.750o-3.375',
 'binding_energies/c-4.750o-2.750',
 'binding_energies/c-4.750o-2.125',
 'binding_energies/c-4.750o-1.500',
 'binding_energies/c-4.062o-6.500',
 'binding_energies/c-4.062o-5.875',
 'binding_energies/c-4.062o-5.250',
 'binding_energies/c-4.062o-4.625',
 'binding_energies/c-4.062o-4.000',
 'binding_energies/c-4.062o-3.375',
 'binding_energies/c-4.062o-2.750',
 'binding_energies/c-4.062o-2.125',
 'binding_energies/c-4.062o-1.500',
 'binding_energies/c-3.375o-6.500',
 'binding_energies/c-3.375o-5.875',
 'binding_energies/c-3.375o-5.250',
 'binding_energies/c-3.375o-4.625',
 'binding_energies/c-3.375o-4.000',
 'binding_energies/c-3.375o-3.375',
 'binding_energies/c-3.375o-2.750',
 'binding_energies/c-3.375o-2.125',
 'binding_energies/c-3.375o-1.500',
 'binding_energies/c-2.688o-6.500',
 'binding_energies/c-2.688o-5.875',
 'binding_energies/c-2.688o-5.250',
 'binding_energies/c-2.688o-4.625',
 'binding_energies/c-2.688o-4.000',
 'binding_energies/c-2.688o-3.375',
 'binding_energies/c-2.688o-2.750',
 'binding_energies/c-2.688o-2.125',
 'binding_energies/c-2.688o-1.500',
 'binding_energies/c-2.000o-6.500',
 'binding_energies/c-2.000o-5.875',
 'binding_energies/c-2.000o-5.250',
 'binding_energies/c-2.000o-4.625',
 'binding_energies/c-2.000o-4.000',
 'binding_energies/c-2.000o-3.375',
 'binding_energies/c-2.000o-2.750',
 'binding_energies/c-2.000o-2.125',
 'binding_energies/c-2.000o-1.500']
In [23]:
# Now run the simulations using a pool.
# Don't run this cell unless you have a while to wait!! ###

#pool = multiprocessing.Pool()
#result = pool.map(run_simulation, experiments)

# Instead, upload the input files to a cluster and run `start_all.sh`
# The script `stage_all.sh` will stage (add) the results in git
# so you can commit and get them back to analyze with the subsequent cells:
In [24]:
base_directory = 'binding_energies'
# base_directory = 'binding_energies_local'
In [25]:
def get_data(experiment):
    carbon, oxygen = experiment
    directory(carbon,oxygen)
    try:
        data = get_pandas_data(directory(carbon,oxygen))
    except OSError:
        return None
    rename_columns(data)
    return data
In [26]:
datas = {tuple(e): get_data(e) for e in experiments}
In [27]:
def get_max_co2(experiment):
    data = datas[tuple(experiment)]
    if data is None:
        return np.nan
    return data[['CO2']].max()
highest_co2 = max([float(get_max_co2(e)) for e in experiments])

Fit an exponential¶

For this first attempt to extract "rate" we fit an exponential growth curve to the normalized CO2 concentration profile of each simulation. First we plot all the curves on one plot, then we'll plot each with its fitted exponential.

In [28]:
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))

ax = plt.axes()
def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for {}".format(experiment))
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values / highest_co2
    ax.plot(times, normalized, label=make_label(experiment))
    #normalized.plot.line(ax=ax, label=make_label(experiment))
    #ax.text(times[-10], normalized[-10], make_label(experiment))

#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO$_2$ concentration')
Out[28]:
Text(0,0.5,'Normalized CO$_2$ concentration')
In [36]:
sns.set_palette('Set1')

def my_function(time, rate):
    "Thing we want to fit."
    return 1. - np.exp(-1*time*rate)
import scipy.optimize


def fit_rate(data):
    normalized = data[['CO2']] / highest_co2
    x_data = np.array(normalized.index)
    y_data = normalized.values[:,0] 
    
    try:
        popt, pcov = scipy.optimize.curve_fit(my_function, x_data, y_data)
    except RuntimeError as e:
        print("☝️ couldn't fit: {}".format(e.message))
        return np.nan
        
    optimal_parameters = popt
    parameter_errors = np.sqrt(np.diag(pcov))
    print("Rate: {} +/- {} (1 st. dev.)".format(optimal_parameters[0],parameter_errors[0]))

    plt.plot(x_data, y_data, 'o')
    plt.plot(x_data, my_function(x_data, *optimal_parameters))
    plt.xlabel('Time (s)')
    plt.ylabel('Normalized CO$_2$ concentration')
    
    plt.show()
    fitted_rate = optimal_parameters[0]
    return fitted_rate

rates = []
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        rates.append(np.nan)
        continue
    rate = fit_rate(data)
    rates.append(rate)

rates
    
[-7.5 -6.5]
Rate: 0.0112469801237 +/- 0.000110181792409 (1 st. dev.)
[-7.5   -5.875]
Rate: 0.247137851008 +/- 0.00115798268134 (1 st. dev.)
[-7.5  -5.25]
Rate: 1.52869663215 +/- 0.0299444056018 (1 st. dev.)
[-7.5   -4.625]
Rate: 299.775535286 +/- 1.64748011008 (1 st. dev.)
[-7.5 -4. ]
Rate: 780.141288463 +/- 2.55214302016 (1 st. dev.)
[-7.5   -3.375]
Rate: 217.137879832 +/- 2.14830908309 (1 st. dev.)
[-7.5  -2.75]
Rate: 3.84229495506 +/- 0.0549194105397 (1 st. dev.)
[-7.5   -2.125]
Rate: 1853.76620184 +/- 20.6659419534 (1 st. dev.)
[-7.5 -1.5]
Rate: 418.10221117 +/- 4.7724621394 (1 st. dev.)
[-6.8125 -6.5   ]
Rate: 0.00890853353773 +/- 4.18081719862e-05 (1 st. dev.)
[-6.8125 -5.875 ]
Rate: 0.238274250807 +/- 0.00117995926324 (1 st. dev.)
[-6.8125 -5.25  ]
Rate: 12.0806242877 +/- 0.133872033318 (1 st. dev.)
[-6.8125 -4.625 ]
Rate: 41.2723298577 +/- 0.215388112267 (1 st. dev.)
[-6.8125 -4.    ]
Rate: 653.371541695 +/- 2.26334521421 (1 st. dev.)
[-6.8125 -3.375 ]
Rate: 195.941035119 +/- 1.72071504296 (1 st. dev.)
[-6.8125 -2.75  ]
Rate: 6.75822526441 +/- 0.0859320768548 (1 st. dev.)
[-6.8125 -2.125 ]
Rate: 1914.1792995 +/- 22.353616955 (1 st. dev.)
[-6.8125 -1.5   ]
Rate: 257.005279522 +/- 3.06461275999 (1 st. dev.)
[-6.125 -6.5  ]
Rate: 0.00310362971828 +/- 1.36980764605e-05 (1 st. dev.)
[-6.125 -5.875]
Rate: 0.0946802693081 +/- 0.000410237476374 (1 st. dev.)
[-6.125 -5.25 ]
Rate: 2.40186777577 +/- 0.0138026119953 (1 st. dev.)
[-6.125 -4.625]
Rate: 1.39779618198 +/- 0.00865924465081 (1 st. dev.)
[-6.125 -4.   ]
Rate: 118.456292223 +/- 0.796709893238 (1 st. dev.)
[-6.125 -3.375]
Rate: 135.073683862 +/- 1.10393343231 (1 st. dev.)
[-6.125 -2.75 ]
Rate: 4.26657279182 +/- 0.0475398030145 (1 st. dev.)
[-6.125 -2.125]
Rate: 1.44629420268 +/- 0.032326506173 (1 st. dev.)
[-6.125 -1.5  ]
Rate: 227.533273841 +/- 2.98091600795 (1 st. dev.)
[-5.4375 -6.5   ]
Rate: 4.86954101342e-05 +/- 7.46199096699e-07 (1 st. dev.)
[-5.4375 -5.875 ]
Rate: 0.00054696547621 +/- 1.03452146354e-05 (1 st. dev.)
[-5.4375 -5.25  ]
Rate: 0.00209549637007 +/- 2.57093330904e-05 (1 st. dev.)
[-5.4375 -4.625 ]
Rate: 0.00169944341792 +/- 1.54372308514e-05 (1 st. dev.)
[-5.4375 -4.    ]
Rate: 3.48228588113 +/- 0.0351064766279 (1 st. dev.)
[-5.4375 -3.375 ]
Rate: 70.4485298251 +/- 0.77339844584 (1 st. dev.)
[-5.4375 -2.75  ]
Rate: 4.8974750312 +/- 0.041252081917 (1 st. dev.)
[-5.4375 -2.125 ]
Rate: 1.3317076701 +/- 0.0209304091091 (1 st. dev.)
[-5.4375 -1.5   ]
Rate: 145.832727001 +/- 1.88082625796 (1 st. dev.)
[-4.75 -6.5 ]
Rate: 9.21095412071e-06 +/- 1.68946370726e-07 (1 st. dev.)
[-4.75  -5.875]
Rate: 1.69920759479e-05 +/- 3.17671242268e-07 (1 st. dev.)
[-4.75 -5.25]
Rate: 1.8482281437e-05 +/- 3.40858214226e-07 (1 st. dev.)
[-4.75  -4.625]
Rate: 1.87327238947e-05 +/- 3.13750869581e-07 (1 st. dev.)
[-4.75 -4.  ]
Rate: 0.388709250369 +/- 0.00419100789149 (1 st. dev.)
[-4.75  -3.375]
Rate: 1.06469374948 +/- 0.0164891098538 (1 st. dev.)
[-4.75 -2.75]
Rate: 0.997624625784 +/- 0.0230269569792 (1 st. dev.)
[-4.75  -2.125]
Rate: 3.13272979337 +/- 0.0683272862321 (1 st. dev.)
[-4.75 -1.5 ]
Rate: 5.2292662113 +/- 0.0210300492995 (1 st. dev.)
[-4.0625 -6.5   ]
Rate: 4.50985600599e-07 +/- 1.06042786122e-08 (1 st. dev.)
[-4.0625 -5.875 ]
Rate: 6.15742238523e-07 +/- 1.45453531882e-08 (1 st. dev.)
[-4.0625 -5.25  ]
Rate: 5.56440543448e-07 +/- 1.23373607929e-08 (1 st. dev.)
[-4.0625 -4.625 ]
Rate: 3.57894293521e-07 +/- 7.55322705229e-09 (1 st. dev.)
[-4.0625 -4.    ]
Rate: 0.233948147032 +/- 0.00111697842629 (1 st. dev.)
[-4.0625 -3.375 ]
Rate: 0.304917480641 +/- 0.00266844148853 (1 st. dev.)
[-4.0625 -2.75  ]
Rate: 0.89127234745 +/- 0.0134413830257 (1 st. dev.)
[-4.0625 -2.125 ]
Rate: 0.305041908696 +/- 0.00382839428718 (1 st. dev.)
[-4.0625 -1.5   ]
Rate: 0.467292043697 +/- 0.00158869489931 (1 st. dev.)
[-3.375 -6.5  ]
Rate: -5.4552238109e-11 +/- inf (1 st. dev.)
[-3.375 -5.875]
Rate: -6.55855553904e-11 +/- inf (1 st. dev.)
[-3.375 -5.25 ]
Rate: -1.26892338488e-10 +/- inf (1 st. dev.)
[-3.375 -4.625]
Rate: -5.41624793145e-11 +/- inf (1 st. dev.)
[-3.375 -4.   ]
Rate: 0.0378554237682 +/- 0.000241607269491 (1 st. dev.)
[-3.375 -3.375]
Rate: 0.579233561548 +/- 0.00226727144247 (1 st. dev.)
[-3.375 -2.75 ]
Rate: 0.600831710121 +/- 0.00346907440536 (1 st. dev.)
[-3.375 -2.125]
Rate: 0.0487209597023 +/- 0.000776848920523 (1 st. dev.)
[-3.375 -1.5  ]
Rate: 0.0373303429542 +/- 0.000121296882531 (1 st. dev.)
[-2.6875 -6.5   ]
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:17: DeprecationWarning: BaseException.message has been deprecated as of Python 2.6
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -5.875 ]
Rate: -1.34770635682e-10 +/- inf (1 st. dev.)
[-2.6875 -5.25  ]
Rate: -2.18187629145e-10 +/- inf (1 st. dev.)
[-2.6875 -4.625 ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -4.    ]
Rate: 0.0145652281822 +/- 6.94430269106e-05 (1 st. dev.)
[-2.6875 -3.375 ]
Rate: 0.329646347808 +/- 0.000360054771534 (1 st. dev.)
[-2.6875 -2.75  ]
Rate: 0.1380947181 +/- 0.000669273531142 (1 st. dev.)
[-2.6875 -2.125 ]
Rate: 0.0110957978742 +/- 0.000239679337853 (1 st. dev.)
[-2.6875 -1.5   ]
Rate: 0.00177576223546 +/- 3.1914747708e-06 (1 st. dev.)
[-2.  -6.5]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.    -5.875]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.   -5.25]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.    -4.625]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2. -4.]
Rate: 0.00235993895126 +/- 1.87566477944e-05 (1 st. dev.)
[-2.    -3.375]
Rate: 0.0244298212732 +/- 3.59759214908e-06 (1 st. dev.)
[-2.   -2.75]
Rate: 0.00151714848026 +/- 9.80487090429e-06 (1 st. dev.)
[-2.    -2.125]
Rate: -2.01387761091e-10 +/- inf (1 st. dev.)
[-2.  -1.5]
Rate: 0.00147317498221 +/- 5.46471742805e-06 (1 st. dev.)
Out[36]:
[0.01124698012372417,
 0.24713785100783817,
 1.5286966321475626,
 299.7755352861337,
 780.1412884625455,
 217.13787983248469,
 3.8422949550577257,
 1853.7662018431663,
 418.1022111704884,
 0.0089085335377349,
 0.23827425080661607,
 12.08062428767228,
 41.27232985771965,
 653.371541694909,
 195.94103511916572,
 6.758225264408982,
 1914.1792994956636,
 257.00527952227026,
 0.0031036297182837,
 0.09468026930806342,
 2.4018677757726885,
 1.3977961819781304,
 118.45629222269837,
 135.07368386180656,
 4.266572791815656,
 1.4462942026759176,
 227.53327384079557,
 4.869541013423381e-05,
 0.0005469654762103199,
 0.0020954963700747244,
 0.0016994434179238584,
 3.4822858811337922,
 70.44852982505572,
 4.8974750311955955,
 1.3317076700954096,
 145.83272700133182,
 9.210954120706608e-06,
 1.6992075947884442e-05,
 1.8482281437005835e-05,
 1.8732723894718023e-05,
 0.3887092503693823,
 1.0646937494834212,
 0.9976246257835009,
 3.132729793367422,
 5.229266211298521,
 4.509856005994856e-07,
 6.157422385234128e-07,
 5.564405434483252e-07,
 3.57894293521307e-07,
 0.23394814703162545,
 0.3049174806413752,
 0.8912723474501224,
 0.30504190869619363,
 0.46729204369717586,
 -5.455223810895606e-11,
 -6.558555539044971e-11,
 -1.2689233848757818e-10,
 -5.4162479314506367e-11,
 0.037855423768156066,
 0.5792335615482348,
 0.6008317101206624,
 0.04872095970227352,
 0.03733034295422652,
 nan,
 -1.3477063568169707e-10,
 -2.1818762914525485e-10,
 nan,
 0.014565228182151937,
 0.3296463478080672,
 0.1380947181001194,
 0.011095797874180795,
 0.0017757622354604488,
 nan,
 nan,
 nan,
 nan,
 0.0023599389512601442,
 0.024429821273206503,
 0.001517148480264792,
 -2.0138776109080676e-10,
 0.0014731749822127736]
In [37]:
rates = np.array(rates)
fixed_rates = rates * (rates>0) + (1e-9 * (rates<0))
log_rates = np.log(fixed_rates)
log_rates
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in greater
  
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in less
  
Out[37]:
array([-4.48765562e+00, -1.39780900e+00,  4.24415498e-01,  5.70303398e+00,
        6.65947504e+00,  5.38053254e+00,  1.34606983e+00,  7.52497463e+00,
        6.03572593e+00, -4.72074564e+00, -1.43433295e+00,  2.49160287e+00,
        3.72019230e+00,  6.48214594e+00,  5.27781377e+00,  1.91076032e+00,
        7.55704425e+00,  5.54909663e+00, -5.77518298e+00, -2.35724965e+00,
        8.76246675e-01,  3.34896841e-01,  4.77454405e+00,  4.90582044e+00,
        1.45081088e+00,  3.69004563e-01,  5.42729649e+00, -9.92992578e+00,
       -7.51112487e+00, -6.16796482e+00, -6.37745448e+00,  1.24768894e+00,
        4.25488237e+00,  1.58871977e+00,  2.86462081e-01,  4.98246026e+00,
       -1.15951171e+01, -1.09827634e+01, -1.08986980e+01, -1.08852386e+01,
       -9.44923643e-01,  6.26871986e-02, -2.37819989e-03,  1.14190476e+00,
        1.65427096e+00, -1.46118304e+01, -1.43004374e+01, -1.44017055e+01,
       -1.48430282e+01, -1.45265578e+00, -1.18771409e+00, -1.15105233e-01,
       -1.18730611e+00, -7.60800856e-01, -2.07232658e+01, -2.07232658e+01,
       -2.07232658e+01, -2.07232658e+01, -3.27398101e+00, -5.46049495e-01,
       -5.09440400e-01, -3.02164596e+00, -3.28794880e+00,             nan,
       -2.07232658e+01, -2.07232658e+01,             nan, -4.22911822e+00,
       -1.10973487e+00, -1.97981547e+00, -4.50118881e+00, -6.33352552e+00,
                   nan,             nan,             nan,             nan,
       -6.04911953e+00, -3.71195071e+00, -6.49092271e+00, -2.07232658e+01,
       -6.52033536e+00])
In [38]:
rate_grid = np.reshape(log_rates, (grid_size,grid_size))
In [39]:
ax = sns.heatmap(rate_grid)
In [40]:
extent = carbon_range + oxygen_range
extent
Out[40]:
(-7.5, -2, -6.5, -1.5)
In [41]:
# Because the center of a corner pixel is in fact the corner of the grid
# we want to stretch the image a little
c_step = mesh[0,1,0]-mesh[0,0,0]
o_step = mesh[1,0,1]-mesh[1,0,0]
carbon_range2 = (carbon_range[0]-c_step/2, carbon_range[1]+c_step/2)
oxygen_range2 = (oxygen_range[0]-c_step/2, oxygen_range[1]+c_step/2)
extent2 = carbon_range2 + oxygen_range2
extent2
Out[41]:
(-7.84375, -1.65625, -6.84375, -1.15625)
In [42]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
plt.plot(-5.997, -4.485, 'ok')
plt.text(-5.997, -4.485, 'Ni(111)')
Out[42]:
Text(-5.997,-4.485,'Ni(111)')
In [43]:
# Binding energies extracted from:
u"""
Medford, A. J.; Lausche, A. C.; Abild-Pedersen, F.;
Temel, B.; Schjødt, N. C.; Nørskov, J. K.; Studt, F. 
Activity and Selectivity Trends in 
Synthesis Gas Conversion to Higher Alcohols. 
Topics in Catalysis 2014, 57 (1-4), 135–142 
DOI: 10.1007/s11244-013-0169-0.
"""
medford_energies = { # Carbon, then Oxygen
'Ru': ( 0.010349288486416697, -2.8153856448231256),
'Rh': ( 0.16558861578266493, -2.546620868091181),
'Ni': ( 0.3001293661060802, -2.5881741535441853),
'Ir': ( 0.36222509702457995, -2.826185484230718),
'Pd': ( 0.28460543337645516, -1.207119596734621),
'Pt': ( 0.8796895213454077, -1.445820136503547),
'Cu': ( 2.323415265200518, -1.7218249542757729),
'Ag': ( 3.855109961190168, -0.8341504215550701),
'Au': ( 3.5601552393272975, -0.10963108355266138),
}
In [44]:
# Shift Medford's energies so that Ni matches Wayne Blaylock's Ni
blaylock_ni = np.array([-5.997, -4.485])
old_ni = np.array(medford_energies['Ni'])
shifted_energies = {metal: tuple(blaylock_ni + np.array(E)-old_ni) for metal,E in medford_energies.items()}
shifted_energies
Out[44]:
{'Ag': (-2.442019404915912, -2.730976268010885),
 'Au': (-2.7369741267787826, -2.006456930008476),
 'Cu': (-3.973714100905562, -3.618650800731588),
 'Ir': (-5.9349042690815, -4.723011330686532),
 'Ni': (-5.997, -4.484999999999999),
 'Pd': (-6.012523932729625, -3.1039454431904363),
 'Pt': (-5.417439844760672, -3.3426459829593616),
 'Rh': (-6.131540750323415, -4.443446714546996),
 'Ru': (-6.286780077619664, -4.712211491278941)}
In [45]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in shifted_energies.iteritems():
    plt.plot(coords[0], coords[1], 'ok')
    plt.text(coords[0], coords[1], metal)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
Out[45]:
(-6.5, -1.5)
In [46]:
# Binding energies for close packed surfaces from
"""
Abild-Pedersen, F.; Greeley, J.; Studt, F.; Rossmeisl, J.;
Munter, T. R.; Moses, P. G.; Skúlason, E.; Bligaard, T.; Norskov, J. K. 
Scaling Properties of Adsorption Energies for Hydrogen-Containing 
Molecules on Transition-Metal Surfaces. 
Phys. Rev. Lett. 2007, 99 (1), 016105 
DOI: 10.1103/PhysRevLett.99.016105.
"""
abildpedersen_energies = { # Carbon, then Oxygen
'Ru': ( -6.397727272727272, -5.104763568600047),
'Rh': ( -6.5681818181818175, -4.609771721406942),
'Ni': ( -6.045454545454545, -4.711681807593758),
'Ir': ( -6.613636363636363, -5.94916142557652),
'Pd': ( -6, -3.517877940833916),
'Pt': ( -6.363636363636363, -3.481481481481482),
'Cu': ( -4.159090909090907, -3.85272536687631),
'Ag': ( -2.9545454545454533, -2.9282552993244817),
'Au': ( -3.7499999999999973, -2.302236198462614),
}
In [47]:
abildpedersen_energies['Pt'][0] - abildpedersen_energies['Ni'][0]
Out[47]:
-0.31818181818181834
In [48]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[48]:
Text(0,0.5,'$\\Delta E^O$ (eV)')

Alternative rate definition¶

For this, we plot the CO2 concentration profiles on a log-log plot, and see the characteristic time as the time at which it crosses the diagonal, i.e. when does $\left( \log_{10}(time) + \log_{10}(\frac{CO_2}{CO_2max}) \right) \ge -10$

In [49]:
"""
Plot everything on a log-log plot
"""
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))
ax = plt.axes()

def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    ax.loglog(times, normalized, label=make_label(experiment))
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
    except IndexError:
        i = -1
    print i, times[i], normalized[i]
    ax.text(times[i], normalized[i], make_label(experiment), rotation=45, ha='center', va='center', fontsize=12)
    
plt.plot([1e-10, 1],[1,1e-10], 'g:')
plt.ylim(1e-10, 1)
plt.xlim(1e-10, 1)
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO2 concentration')
[-7.5 -6.5]
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:23: RuntimeWarning: divide by zero encountered in log10
5546 0.0007183151766476987 1.3926015214032478e-07
[-7.5   -5.875]
4484 2.361167231396705e-05 4.258998833501526e-06
[-7.5  -5.25]
3654 1.6180894144711328e-06 6.264532217051015e-05
[-7.5   -4.625]
3122 1.234189667516712e-06 8.171969781501936e-05
[-7.5 -4. ]
2997 9.750094067547085e-07 0.00010358968343684949
[-7.5   -3.375]
2396 2.961810200801705e-07 0.0003399545734940184
[-7.5  -2.75]
2192 2.4601697082311633e-07 0.00040759368695898393
[-7.5   -2.125]
2758 2.584255024671953e-07 0.00039088982372593805
[-7.5 -1.5]
2730 4.814095903974514e-07 0.00020855548464742454
[-6.8125 -6.5   ]
5207 0.0001974050378455577 5.074959130470027e-07
[-6.8125 -5.875 ]
2970 6.708525102194166e-06 1.4945326031933518e-05
[-6.8125 -5.25  ]
3702 1.9067072157434727e-06 5.285077680775567e-05
[-6.8125 -4.625 ]
3077 2.1391226704174787e-06 4.679404028431258e-05
[-6.8125 -4.    ]
3019 1.192328488927855e-06 8.397665276765415e-05
[-6.8125 -3.375 ]
2604 3.827597685212609e-07 0.0002631705630735117
[-6.8125 -2.75  ]
2003 3.6033935562044284e-07 0.0002787429948988961
[-6.8125 -2.125 ]
2727 3.8056844843973033e-07 0.00026305052512937355
[-6.8125 -1.5   ]
2480 1.3559610986055364e-06 7.43509386574425e-05
[-6.125 -6.5  ]
5217 0.0001254108181920156 7.995409972673579e-07
[-6.125 -5.875]
2838 7.024459996194612e-06 1.4310760084413504e-05
[-6.125 -5.25 ]
3671 5.264207904592093e-06 1.9009161733008402e-05
[-6.125 -4.625]
3602 6.030924843951105e-06 1.666136694638468e-05
[-6.125 -4.   ]
3211 1.9461683273834993e-06 5.149531303727627e-05
[-6.125 -3.375]
2953 6.477543110649682e-07 0.00015508520468369313
[-6.125 -2.75 ]
2325 6.38356184444001e-07 0.00015698837579852698
[-6.125 -2.125]
2529 7.905320486190664e-07 0.0001272044098564433
[-6.125 -1.5  ]
2517 4.769481845257964e-06 2.1009400615745118e-05
[-5.4375 -6.5   ]
4885 0.0001804710079705069 5.556256362329106e-07
[-5.4375 -5.875 ]
2817 1.9548113927223676e-05 5.1206633991158045e-06
[-5.4375 -5.25  ]
3480 2.0628488481728774e-05 4.855892315244647e-06
[-5.4375 -4.625 ]
3662 2.492839986183896e-05 4.01461743967697e-06
[-5.4375 -4.    ]
3073 3.51989168008061e-06 2.8588589807384764e-05
[-5.4375 -3.375 ]
2537 1.138006132023428e-06 8.902019245981242e-05
[-5.4375 -2.75  ]
3072 1.1789193249813288e-06 8.493341148492004e-05
[-5.4375 -2.125 ]
2299 2.529547825931426e-06 3.9904708881754996e-05
[-5.4375 -1.5   ]
2441 1.1586913240120633e-05 8.840470551852554e-06
[-4.75 -6.5 ]
3912 0.00019778430708302343 5.0836952023907e-07
[-4.75  -5.875]
3976 9.832988480042397e-05 1.0222874952334984e-06
[-4.75 -5.25]
3876 0.00012435848709852404 8.125550941902441e-07
[-4.75  -4.625]
4301 0.0001431510856756366 6.990032434397042e-07
[-4.75 -4.  ]
3303 8.995522848481144e-06 1.111838720051659e-05
[-4.75  -3.375]
3274 2.7387210974962524e-06 3.6700334412714816e-05
[-4.75 -2.75]
3154 3.267671285798288e-06 3.0660003627709584e-05
[-4.75  -2.125]
2949 4.8986542783634785e-06 2.0470826838209302e-05
[-4.75 -1.5 ]
2458 2.0826211490285144e-05 4.8141003645429035e-06
[-4.0625 -6.5   ]
4831 0.001156016639277569 8.657616210992914e-08
[-4.0625 -5.875 ]
4665 0.0011185126715017207 8.94285240547611e-08
[-4.0625 -5.25  ]
4630 0.0015421697484169874 6.504143876569925e-08
[-4.0625 -4.625 ]
4770 0.0019194018929889747 5.223061158032518e-08
[-4.0625 -4.    ]
3851 4.920191346933894e-05 2.037549867938948e-06
[-4.0625 -3.375 ]
4026 1.2114270907537852e-05 8.277717523060943e-06
[-4.0625 -2.75  ]
2822 7.723911220345671e-06 1.3088912291108247e-05
[-4.0625 -2.125 ]
2953 7.156570359414204e-06 1.4108498798742125e-05
[-4.0625 -1.5   ]
2851 4.4822227626979134e-05 2.2537164010026197e-06
[-3.375 -6.5  ]
9371 0.047291487577128914 2.1162854170974503e-09
[-3.375 -5.875]
9000 0.049622296539207265 2.015710326696662e-09
[-3.375 -5.25 ]
8727 0.07153330339382889 1.3997309260323104e-09
[-3.375 -4.625]
9476 0.09091354192233976 1.10078906423604e-09
[-3.375 -4.   ]
4304 0.00021901104263880904 4.5714046994526466e-07
[-3.375 -3.375]
3572 4.654369748508949e-05 2.1514506025325027e-06
[-3.375 -2.75 ]
2612 8.919592132129486e-06 1.1305983266344469e-05
[-3.375 -2.125]
3078 1.0660219084454449e-05 9.41455834569704e-06
[-3.375 -1.5  ]
3054 8.53246970804432e-05 1.177541486386267e-06
[-2.6875 -6.5   ]
-1 1.0002090207594585 8.129439687658071e-11
[-2.6875 -5.875 ]
-1 1.0052924298075945 6.12937153063889e-11
[-2.6875 -5.25  ]
-1 1.0017501690656594 5.2702605730034765e-11
[-2.6875 -4.625 ]
-1 1.001632125574509 5.306480730411955e-11
[-2.6875 -4.    ]
3735 0.00029361289585654824 3.4164038560719975e-07
[-2.6875 -3.375 ]
3589 5.83483922600539e-05 1.7151178823411393e-06
[-2.6875 -2.75  ]
2796 1.0136388609746027e-05 9.945738612254589e-06
[-2.6875 -2.125 ]
3331 1.606416521078016e-05 6.336280893421392e-06
[-2.6875 -1.5   ]
3403 0.00017083387563832702 5.874617529670723e-07
[-2.  -6.5]
-1 1.0010451956952058 2.6160937641195612e-12
[-2.    -5.875]
-1 1.002243924015987 2.335251611307557e-12
[-2.   -5.25]
-1 1.0002529187718827 1.0974939528057448e-12
[-2.    -4.625]
-1 1.0005657360240474 1.079009295767126e-12
[-2. -4.]
4512 0.0011587577332387688 8.654645410130304e-08
[-2.    -3.375]
3719 0.0003816737149926224 2.668950601577692e-07
[-2.   -2.75]
3084 7.840789179176985e-05 1.2798261999411759e-06
[-2.    -2.125]
-1 1.0028774685444888 1.7906201939755737e-15
[-2.  -1.5]
3220 0.00030722984873127846 3.290227348978664e-07
Out[49]:
Text(0,0.5,'Normalized CO2 concentration')
In [50]:
"""
An alternative way of defining rates.
"""
new_rates = []
for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for experiment {}".format(experiment))
        new_rates.append(np.nan)
        continue
        
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
        time = times[i]
    except IndexError:
        time = 1.0
    new_rates.append(1./time)
new_log_rates = np.log(np.array(new_rates))
print new_log_rates
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: divide by zero encountered in log10
  app.launch_new_instance()
[ 7.23860212 10.65376938 13.33426448 13.60509594 13.84081872 15.03229502
 15.21786532 15.16865838 14.54654739  8.53025291 11.91213144 13.17013277
 13.05511478 13.63960245 14.77585828 14.8362196  14.78159979 13.51100006
  8.98391566 11.86611221 12.15457987 12.01861019 13.14964808 14.24975436
 14.26436943 14.05055964 12.25327289  8.61994041 10.84263175 10.7888375
 10.59950285 12.55708034 13.68623283 13.65091237 12.88747    11.36563427
  8.52833348  9.22718256  8.99234214  8.85160994 11.61878357 12.8080195
 12.63143297 12.22655003 10.7792982   6.76277511  6.79575545  6.47456493
  6.25574166  9.91957804 11.32112639 11.77118969 11.84747969 10.01280639
  3.05142497  3.00331502  2.63759216  2.39784631  8.42638841  9.97511896
 11.62726034 11.44899159  9.36904661  0.          0.          0.
  0.          8.13324834  9.74907875 11.49937878 11.03891953  8.67481896
  0.          0.          0.          0.          6.76040677  7.87094446
  9.45358598  0.          8.0879144 ]
In [51]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[51]:
Text(0,0.5,'$\\Delta E^O$ (eV)')
In [52]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='spline16', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[52]:
Text(0,0.5,'$\\Delta E^O$ (eV)')

Count the reactions¶

In [53]:
def get_reaction_count(experiment):
    d = directory(*experiment)
    f = os.path.join(d,'chemkin','chem_annotated.inp')
    with open(f) as chemkin:
        r = re.compile('\! Reaction index: Chemkin \#(\d+)')
        for line in chemkin:
            m = r.match(line)
            if m:
                count = int(m.group(1))
    return count
    
get_reaction_count(experiment)
Out[53]:
373
In [54]:
i=35
print experiments[i]
print get_reaction_count(experiments[i])
[-5.4375 -1.5   ]
536
In [55]:
reaction_counts = map(get_reaction_count, experiments)
reaction_counts_grid = np.log10(np.reshape(reaction_counts, (grid_size,grid_size)))

plt.imshow(reaction_counts_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    color='k'
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range2)
plt.ylim(oxygen_range2)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')

for e,n in zip(experiments,reaction_counts):
    plt.text(e[0],e[1],n,color='w',ha='center', va='center')
#plt.colorbar()
In [56]:
# A linear one, just to check it looks the same
reaction_counts_grid = np.reshape(reaction_counts, (grid_size,grid_size))
ax = sns.heatmap(reaction_counts_grid.T[::-1,:], annot=True, fmt='d', square=True)
In [ ]: